#define CONFIG_USE_DOUBLE 1
#if CONFIG_USE_DOUBLE

#if defined(cl_khr_fp64)  // Khronos extension available?
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#elif defined(cl_amd_fp64)  // AMD extension available?
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#endif

// double
typedef double real_t;
// mass definitions
#define C12_MASS 12.0000000
#define C13_MASS 13.0033548378
#define C12_INT_MASS 12
#define C13_INT_MASS 13
#define H1_MASS 1.0078250321
#define H2_MASS 2.0141017780
#define H1_INT_MASS 1
#define H2_INT_MASS 2
#define N14_MASS 14.0030740052
#define N15_MASS 15.0001088984
#define N14_INT_MASS 14
#define N15_INT_MASS 15
#define O16_MASS 15.9949146221
#define O17_MASS 16.99913150
#define O18_MASS 17.9991604
#define O16_INT_MASS 16
#define O17_INT_MASS 17
#define O18_INT_MASS 18
#define P31_MASS 30.97376151
#define P31_INT_MASS 31
#define S32_MASS 31.97207069
#define S33_MASS 32.97145850
#define S34_MASS 33.96786683
#define S36_MASS 35.96708088
#define S32_INT_MASS 32
#define S33_INT_MASS 33
#define S34_INT_MASS 34
#define S36_INT_MASS 36
#define ELECTRON_MASS 0.00054857990946
#define PROTON_MASS   1.007276466812
#define PI      3.14159265358979323846
#define TWOPI   6.28318530717958647
#define HALFPI  1.57079632679489666

#else

// float
typedef float real_t;
typedef float2 real2_t;
// mass definitions
#define C12_MASS 12.0000000f
#define C13_MASS 13.0033548378f
#define C12_INT_MASS 12
#define C13_INT_MASS 13
#define H1_MASS 1.0078250321f
#define H2_MASS 2.0141017780f
#define H1_INT_MASS 1
#define H2_INT_MASS 2
#define N14_MASS 14.0030740052f
#define N15_MASS 15.0001088984f
#define N14_INT_MASS 14
#define N15_INT_MASS 15
#define O16_MASS 15.9949146221f
#define O17_MASS 16.99913150f
#define O18_MASS 17.9991604f
#define O16_INT_MASS 16
#define O17_INT_MASS 17
#define O18_INT_MASS 18
#define P31_MASS 30.97376151f
#define P31_INT_MASS 31
#define S32_MASS 31.97207069f
#define S33_MASS 32.97145850f
#define S34_MASS 33.96786683f
#define S36_MASS 35.96708088f
#define S32_INT_MASS 32
#define S33_INT_MASS 33
#define S34_INT_MASS 34
#define S36_INT_MASS 36
#define ELECTRON_MASS 0.00054857990946f
#define PROTON_MASS   1.007276466812f
#define PI      3.14159265358f
#define TWOPI   6.28318530717f
#define HALFPI  1.57079632679f

#endif

struct clcomplex {
    real_t real;
    real_t imaginary;
};

real_t pow_pow(real_t a, int b) {
    real_t result = a;
    for (int i = 1; i < b; i++) {
        result *= a;
    }
    return result;
}

struct input_data {
    int id;
    int mode;
    real_t C12;
    real_t C13;
    real_t H1;
    real_t H2;
    real_t N14;
    real_t N15;
    real_t O16;
    real_t O17;
    real_t O18;
    real_t P31;
    real_t S32;
    real_t S33;
    real_t S34;
    real_t S36;
    int numC;
    int numH;
    int numN;
    int numO;
    int numP;
    int numS;
    float accuracy;
};

struct gpu_input_data {
    real_t molecularWeight;
    real_t minMolecularWeight;
    int maxIntIsoShiftAsInt;
    real_t integerMolecularWeight;
    int integerMolecularWeightAsInt;
    int minIntegerMolecularWeightAsInt;
    int maxIntegerMolecularWeightAsInt;
    real_t molecularVariance;
    real_t integerMolecularVariance;
    long massRange;
    long numPoints;
    int resolution;
    int bitSize;
};

__kernel void calculateFrequency(__constant struct gpu_input_data* gpuInput, __constant struct input_data* input, __global struct clcomplex* frequencyData) {
    real_t real, imag, freq, theta, r;

    ulong i2 = get_global_id(0);
    // Frequency Domain calculation for every Point in the output array!
    // Normally one would just compute the Mass Range and then do the permutation of
    // data for the fft (reverse the bit order of the input array), but since it will run on the
    // gpu, we will compute the data of the bit reversed array field and put it in the corect one.
    // this will make the overall computation faster since no additional memory access is needed to swap values
    // reverse long completely, then shift to the right
    ulong temp = i2;                                    // temporary variable needed to do the shift
    ulong i = i2;                                       // i will be reversed bits of temp; 
                                                        // reverse order of every bit
    int s = sizeof (i) * CHAR_BIT - 1;                  // extra shift needed at end

    for (temp >>= 1; temp; temp >>= 1)                  //as long as temp is not zero
    {
        i <<= 1;                                        // shift bits of i to the left
        i |= temp & 1;                                  // if lsb of temp is 1, add it to the right of i
        s--;
    }
    i <<= s;                                            // shift when v's highest bits are zero
    i >>= sizeof (i) * CHAR_BIT - gpuInput->bitSize;    // shift it to the right again,
                                                        // because the variable can hold more bits than needed

    // The upper and lower half of the Points are handled differently
    if (i < gpuInput->numPoints / 2) {
        freq = (real_t) i / gpuInput->massRange;
    } else {
        freq = -(real_t) (gpuInput->numPoints - i) / gpuInput->massRange;
    }


    r = 1;      //radius
    theta = 0;  //angle

    // Do this for every element. Write efficient code for GPU use later on

    // Carbon
    if (input->numC) {
        real = input->C12 * cos(TWOPI * C12_INT_MASS * freq)
                + input->C13 * cos(TWOPI * C13_INT_MASS * freq);
        imag = input->C12 * sin(TWOPI * C12_INT_MASS * freq)
                + input->C13 * sin(TWOPI * C13_INT_MASS * freq);
        r *= pow_pow(hypot(real, imag), input->numC);
        theta = theta + input->numC * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    // Helium
    if (input->numH) {
        real = input->H1 * cos(TWOPI * H1_INT_MASS * freq)
                + input->H2 * cos(TWOPI * H2_INT_MASS * freq);
        imag = input->H1 * sin(TWOPI * H1_INT_MASS * freq)
                + input->H2 * sin(TWOPI * H2_INT_MASS * freq);
        r *= pow_pow(hypot(real, imag), input->numH);
        theta = theta + input->numH * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    // Nitrogen
    if (input->numN) {
        real = input->N14 * cos(TWOPI * N14_INT_MASS * freq)
                + input->N15 * cos(TWOPI * N15_INT_MASS * freq);
        imag = input->N14 * sin(TWOPI * N14_INT_MASS * freq)
                + input->N15 * sin(TWOPI * N15_INT_MASS * freq);
        r *= pow_pow(hypot(real, imag), input->numN);
        theta = theta + input->numN * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    // Oxygen
    if (input->numO) {
        real = input->O16 * cos(TWOPI * O16_INT_MASS * freq)
                + input->O17 * cos(TWOPI * O17_INT_MASS * freq)
                + input->O18 * cos(TWOPI * O18_INT_MASS * freq);
        imag = input->O16 * sin(TWOPI * O16_INT_MASS * freq)
                + input->O17 * sin(TWOPI * O17_INT_MASS * freq)
                + input->O18 * sin(TWOPI * O18_INT_MASS * freq);
        r *= pow_pow(hypot(real, imag), input->numO);
        theta = theta + input->numO * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    // Phosporus
    if (input->numP) {
        real = input->P31 * cos(TWOPI * P31_INT_MASS * freq);
        imag = input->P31 * sin(TWOPI * P31_INT_MASS * freq);
        r *= pow_pow(hypot(real, imag), input->numP);
        theta = theta + input->numP * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    // Sulfur
    if (input->numS) {
        real = input->S32 * cos(TWOPI * S32_INT_MASS * freq)
                + input->S33 * cos(TWOPI * S33_INT_MASS * freq)
                + input->S34 * cos(TWOPI * S34_INT_MASS * freq)
                + input->S36 * cos(TWOPI * S36_INT_MASS * freq);
        imag = input->S32 * sin(TWOPI * S32_INT_MASS * freq)
                + input->S33 * sin(TWOPI * S33_INT_MASS * freq)
                + input->S34 * sin(TWOPI * S34_INT_MASS * freq)
                + input->S36 * sin(TWOPI * S36_INT_MASS * freq);
        r *= pow_pow(hypot(real, imag), input->numS);
        theta = theta + input->numS * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    /* Convert back to real:imag coordinates and store */

    frequencyData[i2].real = r * cos(theta - TWOPI * gpuInput->integerMolecularWeightAsInt * freq); //taken from midas
    frequencyData[i2].imaginary = r * sin(theta - TWOPI * gpuInput->integerMolecularWeightAsInt * freq); //taken from midas
}

// Compute T x DFT-2
// T is the number of threads.
// N = 2*T is the size of input vectors.
// P is the length of input sub-sequences: 1,2,4,...,T.

__kernel void fftRadix2Kernel(__global struct clcomplex* indata, int d) {
    int threadId = get_global_id(0); // global id
    int i = threadId % d + (threadId / d)*(d << 1); // first index
    int j = i + d; // second index
    
    // use "normal" math functions
    real_t e = cospi((real_t) i / (real_t) d);
    real_t f = -sinpi((real_t) i / (real_t) d);
    // use optimized native math
    //real_t e = native_cos(M_PI*(real_t) i / (real_t) d);
    //real_t f = -native_sin(M_PI*(real_t) i / (real_t) d);

    struct clcomplex u0 = indata[i];
    struct clcomplex u1 = indata[j];
    struct clcomplex swap;

    swap.real = (e * u1.real) - (f * u1.imaginary);
    swap.imaginary = (f * u1.real) + (e * u1.imaginary);
    u1.real = u0.real - swap.real;
    u1.imaginary = u0.imaginary - swap.imaginary;
    u0.real += swap.real;
    u0.imaginary += swap.imaginary;

    indata[i] = u0;
    indata[j] = u1;
}

__kernel void calculateFrequencyand32Items(__constant struct gpu_input_data* gpuInput, __constant struct input_data* input, __global struct clcomplex* frequencyData) {
        real_t real, imag, freq_0, freq_1, theta_0,theta_1, r_0, r_1;
        ulong i2 = get_global_id(0)*2;
        ulong temp = i2;                                    // temporary variable needed to do the shift
        ulong i_0=i2;
        ulong i_1=i2+1;
        
        int s = sizeof (i2) * CHAR_BIT - 1;                  // extra shift needed at end

        for (temp >>= 1; temp; temp >>= 1)                  //as long as temp is not zero
        {
            i_0 <<= 1;                                        // shift bits of i to the left
            i_1 <<= 1;                                        // shift bits of i to the left
            i_0 |= temp & 1;                                  // if lsb of temp is 1, add it to the right of i
            i_1 |= temp & 1;                                  // if lsb of temp is 1, add it to the right of i
            s--;
        }                                         // shift when temp's highest bits are zero
        i_0 <<= s;   
        i_1 <<= s;   
        i_0 >>= sizeof (i2) * CHAR_BIT - gpuInput->bitSize;    // shift it to the right again,
        i_1 >>= sizeof (i2) * CHAR_BIT - gpuInput->bitSize;    // shift it to the right again,
        
        
            // The upper and lower half of the Points are handled differently
        if (i_0 < gpuInput->numPoints / 2) {
            freq_0 = (real_t) i_0 / gpuInput->massRange;
        } else {
            freq_0 = -(real_t) (gpuInput->numPoints - i_0) / gpuInput->massRange;
        }
        
        if (i_1 < gpuInput->numPoints / 2) {
            freq_1 = (real_t) i_1 / gpuInput->massRange;
        } else {
            freq_1 = -(real_t) (gpuInput->numPoints - i_1) / gpuInput->massRange;
        }
        
    r_0 = 1;      //radius
    r_1 = 1;      //radius
    theta_0 = 0;  //angle
    theta_1 = 0;  //angle

    // Do this for every element. Write efficient code for GPU use later on

    // Carbon
    if (input->numC) {
        real = input->C12 * cos(TWOPI * C12_INT_MASS * freq_0)
                + input->C13 * cos(TWOPI * C13_INT_MASS * freq_0);
        imag = input->C12 * sin(TWOPI * C12_INT_MASS * freq_0)
                + input->C13 * sin(TWOPI * C13_INT_MASS * freq_0);
        r_0 *= pow_pow(hypot(real, imag), input->numC);
        theta_0 += input->numC * atan2(imag, real); 
        
        real = input->C12 * cos(TWOPI * C12_INT_MASS * freq_1)
                + input->C13 * cos(TWOPI * C13_INT_MASS * freq_1);
        imag = input->C12 * sin(TWOPI * C12_INT_MASS * freq_1)
                + input->C13 * sin(TWOPI * C13_INT_MASS * freq_1);
        r_1 *= pow_pow(hypot(real, imag), input->numC);
        theta_1 += input->numC * atan2(imag, real); //[-pi,pi];	 , taken from midas
    }

    // Helium
    if (input->numH) {
        real = input->H1 * cos(TWOPI * H1_INT_MASS * freq_0)
                + input->H2 * cos(TWOPI * H2_INT_MASS * freq_0);
        imag = input->H1 * sin(TWOPI * H1_INT_MASS * freq_0)
                + input->H2 * sin(TWOPI * H2_INT_MASS * freq_0);
        r_0 *= pow_pow(hypot(real, imag), input->numH);
        theta_0 += input->numH * atan2(imag, real); 
        
        real = input->H1 * cos(TWOPI * H1_INT_MASS * freq_1)
                + input->H2 * cos(TWOPI * H2_INT_MASS * freq_1);
        imag = input->H1 * sin(TWOPI * H1_INT_MASS * freq_1)
                + input->H2 * sin(TWOPI * H2_INT_MASS * freq_1);
        r_1 *= pow_pow(hypot(real, imag), input->numH);
        theta_1 += input->numH * atan2(imag, real); 
    }

    // Nitrogen
    if (input->numN) {
        real = input->N14 * cos(TWOPI * N14_INT_MASS * freq_0)
                + input->N15 * cos(TWOPI * N15_INT_MASS * freq_0);
        imag = input->N14 * sin(TWOPI * N14_INT_MASS * freq_0)
                + input->N15 * sin(TWOPI * N15_INT_MASS * freq_0);
        r_0 *= pow_pow(hypot(real, imag), input->numN);
        theta_0 += input->numN * atan2(imag, real);
        
        real = input->N14 * cos(TWOPI * N14_INT_MASS * freq_1)
                + input->N15 * cos(TWOPI * N15_INT_MASS * freq_1);
        imag = input->N14 * sin(TWOPI * N14_INT_MASS * freq_1)
                + input->N15 * sin(TWOPI * N15_INT_MASS * freq_1);
        r_1 *= pow_pow(hypot(real, imag), input->numN);
        theta_1 += input->numN * atan2(imag, real);
    }

    // Oxygen
    if (input->numO) {
        real = input->O16 * cos(TWOPI * O16_INT_MASS * freq_0)
                + input->O17 * cos(TWOPI * O17_INT_MASS * freq_0)
                + input->O18 * cos(TWOPI * O18_INT_MASS * freq_0);
        imag = input->O16 * sin(TWOPI * O16_INT_MASS * freq_0)
                + input->O17 * sin(TWOPI * O17_INT_MASS * freq_0)
                + input->O18 * sin(TWOPI * O18_INT_MASS * freq_0);
        r_0 *= pow_pow(hypot(real, imag), input->numO);
        theta_0 += input->numO * atan2(imag, real);
        
        real = input->O16 * cos(TWOPI * O16_INT_MASS * freq_1)
                + input->O17 * cos(TWOPI * O17_INT_MASS * freq_1)
                + input->O18 * cos(TWOPI * O18_INT_MASS * freq_1);
        imag = input->O16 * sin(TWOPI * O16_INT_MASS * freq_1)
                + input->O17 * sin(TWOPI * O17_INT_MASS * freq_1)
                + input->O18 * sin(TWOPI * O18_INT_MASS * freq_1);
        r_1 *= pow_pow(hypot(real, imag), input->numO);
        theta_1 += input->numO * atan2(imag, real);
    }

    // Phosporus
    if (input->numP) {
        real = input->P31 * cos(TWOPI * P31_INT_MASS * freq_0);
        imag = input->P31 * sin(TWOPI * P31_INT_MASS * freq_0);
        r_0 *= pow_pow(hypot(real, imag), input->numP);
        theta_0 += input->numP * atan2(imag, real);
        
        real = input->P31 * cos(TWOPI * P31_INT_MASS * freq_1);
        imag = input->P31 * sin(TWOPI * P31_INT_MASS * freq_1);
        r_1 *= pow_pow(hypot(real, imag), input->numP);
        theta_1 += input->numP * atan2(imag, real);
    }

    // Sulfur
    if (input->numS) {
        real = input->S32 * cos(TWOPI * S32_INT_MASS * freq_0)
                + input->S33 * cos(TWOPI * S33_INT_MASS * freq_0)
                + input->S34 * cos(TWOPI * S34_INT_MASS * freq_0)
                + input->S36 * cos(TWOPI * S36_INT_MASS * freq_0);
        imag = input->S32 * sin(TWOPI * S32_INT_MASS * freq_0)
                + input->S33 * sin(TWOPI * S33_INT_MASS * freq_0)
                + input->S34 * sin(TWOPI * S34_INT_MASS * freq_0)
                + input->S36 * sin(TWOPI * S36_INT_MASS * freq_0);
        r_0 *= pow_pow(hypot(real, imag), input->numS);
        theta_0 += input->numS * atan2(imag, real); 
        
        real = input->S32 * cos(TWOPI * S32_INT_MASS * freq_1)
                + input->S33 * cos(TWOPI * S33_INT_MASS * freq_1)
                + input->S34 * cos(TWOPI * S34_INT_MASS * freq_1)
                + input->S36 * cos(TWOPI * S36_INT_MASS * freq_1);
        imag = input->S32 * sin(TWOPI * S32_INT_MASS * freq_1)
                + input->S33 * sin(TWOPI * S33_INT_MASS * freq_1)
                + input->S34 * sin(TWOPI * S34_INT_MASS * freq_1)
                + input->S36 * sin(TWOPI * S36_INT_MASS * freq_1);
        r_1 *= pow_pow(hypot(real, imag), input->numS);
        theta_1 += input->numS * atan2(imag, real); 
    }

    // define variables needed for fft
    struct clcomplex u0;
    struct clcomplex u1;
    struct clcomplex swap;
    
    // safe results in registers
    u0.real = r_0 * cos(theta_0 - TWOPI * gpuInput->integerMolecularWeightAsInt * freq_0); 
    u0.imaginary = r_0 * sin(theta_0 - TWOPI * gpuInput->integerMolecularWeightAsInt * freq_0); 
    u1.real = r_1 * cos(theta_1 - TWOPI * gpuInput->integerMolecularWeightAsInt * freq_1); 
    u1.imaginary = r_1 * sin(theta_1 - TWOPI * gpuInput->integerMolecularWeightAsInt * freq_1); 
    
    // create clcomplex array of size 2*32 (warp size * 2 indizes per thread) = 64 in local memory
    __local struct clcomplex localData[64];
    
    
    for (int d=1; d<64; d<<=1) {
        int threadId = get_local_id(0); // local id
        int i_local = threadId % d + (threadId / d)*(d << 1); // first index in local array
        int j_local = i_local + d; // second index in local array
        
        // use "normal" math functions
        real_t e = cospi((real_t) i_local / (real_t) d);
        real_t f = -sinpi((real_t) i_local / (real_t) d);
        
        if(d!=1) {
           u0 = localData[i_local];
           u1 = localData[j_local];
        }
        
        swap.real = (e * u1.real) - (f * u1.imaginary);
        swap.imaginary = (f * u1.real) + (e * u1.imaginary);
        u1.real = u0.real - swap.real;
        u1.imaginary = u0.imaginary - swap.imaginary;
        u0.real += swap.real;
        u0.imaginary += swap.imaginary;
        
           localData[i_local]=u0;
           localData[j_local]=u1;
        
    }
    
    frequencyData[get_global_id(0)]=localData[get_local_id(0)];
    frequencyData[get_global_id(0)+32]=localData[get_local_id(0)+32];
}